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NONLINEAR  EQUATORIAL  SPREAF  F: 

THE  EFFECT  OF  NEUTRAL  WINDS  AND 
BACKGROUND  PEDERSEN  CONDUCTIVITY 

1.  Introduction 

In  our  previous  studies  of  evolving  equatorial  spread  F  (ESF)  bubbles 
and  plumes  in  the  equatorial  ionosphere  f  Scannapieco  and  Ossakow,  1976; 
Ossakow  et  al. ,  1979;  Zalesak  and  Ossakow.  1980}  ,  the  effects  of  the  neutral 
wind  were  neglected.  Rather,  we  focused  our  attention  on  showing  that  the 
motion  and  structure  of  the  experimentally  measured  ESF  environment  (bottom- 
side  and  topside  spread  F,  bubble  formation  and  evolution)  could  be  ex¬ 
plained  in  terms  of  the  nonlinear  evolution  of  the  gravitationally  driven 
collisional  Rayleigh- Taylor  instability.  Through  the  use  of  numerical  simu¬ 
lation  techniques  we  were  able  to  demonstrate  both  qualitative  and  quantita¬ 
tive  agreement  with  the  observations.  We  wished  to  show  that  a  simple  model 
(i.e. ,  using  only  gravity  and  the  bottomside  background  electron  density 
gradient  as  drivers),  followed  into  the  nonlinear  regime,  could  explain 
observations  that  were  up  to  that  point  inexplicable.  However,  there  are 
some  aspects  of  the  observations  which  we  do  not  see  in  our  previous  simu¬ 
lations. 

First,  there  is  the  tendency  of  ESF  structure  to  drift  eastward  at 
approximately  the  neutral  wind  velocity,  obviously  something  which  could 
not  be  duplicated  in  a  numerical  simulation  which  neglected  neutral  wind 
effects.  Secondly,  there  is  the  curious  tendency  on  the  part  of  radar  back- 
scatter  maps  of  ESF  to  show  plumes  of  backscatter  intensity  which  tilt  east¬ 
ward  with  altitude  at  the  lower  altitudes  and  westward  with  altitude  at  the 
higher  altitudes.  These  structures  were  dubbed  "C's"  and  "fishtails"  by 
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Woodman  and  La  Hoz  [1976]  and  have  also  been  seen  In  the  ALTAIR  backscatter 
maps  of  Tsunoda  [l98l] ,  although  Tsunoda  chooses  not  to  regard  the  eastward- 
tilting  and  westward-tilting  structures  as  part  of  the  same  plume.  Addition¬ 
ally,  we  should  point  out  that  McClure  et  al.,  [1977]  have  observed  the 
westward  drift  of  plasma  bubbles  (bite-outs  or  depletions).  We  propose  here 
a  simple  model  of  the  interaction  of  the  eastward  neutral  wind  at  the  equator 
with  the  equatorial  ionosphere  which  we  believe  explains  both  of  these  ob¬ 
servations.  At  this  juncture,  we  should  point  out  that  Woodman  and  La  Hoz 
Q 1976 ]  ,  Ott  f 1978],  and  Ossakow  and  Chaturvedi  [1978]  hypothesized  that  an 
eastward  neutral  wind  would  produced  westward  drift  of  ESF  bubbles. 

Briefly,  we  find  that  if  the  magnetic  field  line  integrated  Pedersen 
conductivity  has  a  finite  contribution  from  plasma  which  is  not  subject  to 
the  equatorial  F  region  neutral  wl~nd  (e.g.,  plasma  at  higher  latitude  E 
regions),  then  the  vertical  polarization  electric  field  driven  by  the  neutral 
wind  at  the  equator  is  partially  shorted  out  by  this  "background"  E  region 
conductivity,  causing  there  to  be  relative  motion,  or  "slip",  between  the 
plasma  and  the  neutral  wind  at  the  equator.  This  effect  was  first  described 
by  Rishbeth  [ 19 7 1 ]  .  Further  investigation  shows  that  the  degree  of  "slip" 
is  inversely  proportional  to  the  "local"  (i.e. ,  equatorial  F  region)  Pedersen 
conductivity,  causing  there  to  be  a  vertical  shear  in  the  plasma  motion  even 
when  there  is  no  vertical  shear  in  the  neutral  wind.  This  plasma  shear  bends 
any  vertical  structure  about  the  "local"  maximum  in  Pedersen  conductivity, 
giving  rise  to  the  "C's"  and  "fishtails"  seen  on  coherent  backscatter  radar 
maps  (Woodman  and  La  Hoz,  1976;  Tsunoda ,  1981).  Zalesak  et  al.  (1980)  have 
presented  a  preliminary  version  of  this  model. 

In  section  2  we  present  the  geometry  of  the  physical  problem  we  are 
modeling  and  briefly  review  the  relevant  equations  of  motion.  We  also  show 
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that  any  passive  structure  placed  in  the  ambient  equatorial  environment  will, 
in  fact,  be  bent  into  C-shaped  structures.  However,  spread  F  plumes  are 
far  from  passive  structures,  and  it  is  necessary  to  perform  nonlinear  numeric 
cal  simulations  to  prove  the  case  unequivocally.  These  simulations  are 
presented  in  section  3,  where  we  also  show  the  surprising  result  that,  for 
the  case  studied,  it  may  be  the  eastward  wall  of  the  plume,  as  well  as  the 
westward  wall,  which  is  subject  to  secondary  instabilities.  A  stability 
analysis  which  included  only  the  interaction  of  the  neutral  wind  with  the 
plasma  gradients  in  the  bubble  would  conclude  that  it  should  be  only  the 
westward  wall  of  the  plume  which  is  unstable.  Consideration  of  the  self- 
consistent  polarization  electric  field  of  the  bubble  itself,  as  well  as  of 
gravitational  effects  on  the  tilted  structure  can  cause  the  instability  to 
"switch  sides".  In  section  4  we  present  our  conclusions,  and  in  section  5 
we  discuss  briefly  our  plans  for  future  work. 
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2.  Theory 


In  Figure  1,  we  show  the  geometry  of  the  physical  phenomenon  we  are 
attempting  to  model.  The  equatorial  F  region  plasma  responds  to  the  effects 
of  the  earth’s  magnetic  field,  gravity,  collisions  with  the  neutral  atmosphere, 
and  electric  fields.  Since  the  conductivity  along  magnetic  field  lines  is 
extremely  high,  these  electric  fields  can  depend  on  the  dynamics  of  plasma 
far  from  the  equatorial  region,  but  connected  to  the  equatorial  region  by 
magnetic  field  lines.  We  find  that  the  physical  quantity  dominating  the 
evolution  of  the  collisional  Rayleigh-Taylor  instability  is  the  magnetic 
field  line  integrated  Pedersen  conductivity,  and  that  the  primary  contri¬ 
bution  to  that  quantity  comes  from  plasma  in  the  local  region  near  the 
"computational  plane"  shown  in  Figure  1.  This  fact  has  been  the  basis  for 
our  previous  theoretical  and  numerical  studies  of  equatorial  spread  F, 
(Scannapieco  and  Ossakow,  1976;  Ossakow  et  al. ,  1979;  Zalesak  and  Ossakow, 

1980)  and  has  enabled  us  to  study  the  phenomena  of  interest  using  just  a 
single  two-dimensional  computational  plane. 

We  do  not  propose  here  to  analyze  the  problem  in  the  complete  three- 
dimensional  geometry,  but  rather,  as  a  first  step,  to  modify  our  two- 
dimensional  model  to  take  into  account  the  presence  of  other  plasma,  and 
hence  Pedersen  conductivities  and  forces,  in  regions  far  from  the  equatorial 
plane,  but  connected  to  the  equatorial  F  region  plasma  along  magnetic  field 
lines.  For  instance  this  could  be  the  northern  and  southern  hemisphere  E 
region  plasma  shown  in  Figure  1.  This  modification  is  shown  in  Figure  2, 
where  we  show  three  distinct  layers  of  plasma  connected  by  magnetic  field 
lines.  The  center  layer  is  the  same  computational  plane  as  we  have  used  in 
our  previous  work  (Scannapieco  and  Ossakow,  1976;  Ossakow  et  al.,  1979; 
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EQUATORIAL 
SPREAD  F  MODEL 


.  DASHED  LINE  DEPICTS 
E  REGION  PLASMA 


COMPUTATIONAL  PLANE 
'(EQUATORIAL  F  REGION 
PLASMA) 


Fig.  1  —  Diagram  of  the  equatorial  ionosphere  and  of  the  neighboring  regions  which  have 
physical  relevance  to  equatorial  spread  F  (ESF)  processes,  including  the  E  region  plasma 
at  higher  and  lower  latitudes.  These  regions  are  electrically  coupled  to  the  equatorial  F 
region  ionosphere  by  the  high  conductivity  along  magnetic  field  lines.  Plasma  is  actually 
distributed  all  along  these  field  lines,  but  in  this  study  we  shall  make  the  assumption  that 
this  system  can  be  modeled  accurately  by  three  planes  of  plasma  connected  by  straight 
field  lines,  as  shown  in  Figure  2.  One  of  these  three  layers  (layer  2  in  Figure  2)  is  shown 
here  as  the  “computational  plane.” 
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Fig.  2  —  The  “three  layer”  model  of  the  physical  system  depicted  in  Figure  1.  All  plas¬ 
ma  in  the  vicinity  of  the  equatorial  plane  has  been  compressed  into  layer  2,  while  the 
remaining  northern  and  southern  hemisphere  plasma  has  been  compressed  into  layers  1 
and  3  respectively.  Further,  the  magnetic  field  lines  have  been  straightened  so  we  can 
deal  in  cartesian  coordinates  x,  y,  and  z  as  shown  in  the  figure.  The  plasma  in  layers  1 
and  3  is  assumed  to  be  uniform  and  free  of  any  external  driving  force  such  as  a  neutral 
wind.  The  equatorial  layer  2  is  assigned  a  realistic  initial  distribution  of  electron  den¬ 
sity  N0(y),  and  ion-neutral  collision  frequency,  along  with  a  neutral  wind  which  may 
vary  with  altitude,  but  which  is  taken  to  be  uniform  and  eastward,  and  equal  to 
150  m/sec  in  this  study.  In  addition,  gravity  points  in  the  negative  y  direction. 
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Zalesak  and  Ossakow,  1980)  and  represents  the  equatorial  nighttime  F  region 
plasma.  The  upper  and  lower  layers  represent  the  remaining  northern  and 


southern  hemisphere  plasma  respectively,  including  the  E  region  plasma.  The 
problem  is  still  essentially  two-dimensional  in  that  we  do  not  allow  trans¬ 
port  of  ions  between  layers,  nor  do  we  allow  any  physical  quantity  to  vary 
with  z  within  a  layer,  where  z  is  the  direction  along  the  magnetic  field. 

We  do,  however,  allow  electron  currents  to  flow  along  field  lines  between  the 
layers  to  preserve  electrical  neutrality.  Also,  within  the  context  of  this 
model,  we  will  finally  take  the  E  region  layers  to  act  as  a  passive  load, 
i.e.,  we  do  not  allow  for  any  change  in  layers  1  and  3  and  those  layers  are 
assumed  to  remain  uniform.  Thus,  as  a  first  cut  we  are  taking  our  previous 
equatorial  plane  simulations  (Scannapieco  and  Ossakow,  1976;  Ossakow  et  al. , 
1979;  Zalesak  and  Ossakow,  1980)  and  adding  a  passive  E  region  load  to  the 
circuit  to  allow  for  short  circuiting  ef Sects.  Under  the  assumptions  that: 

(i)  the  electric  fields  of  interest  are  electrostatic  and,  hence,  derivable 
from  a  scalar  potential;  and  (ii)  the  conductivity  along  magnetic  field  lines 
is  extremely  large  and,  hence,  the  potential  is  constant  along  a  field  line, 
we  are  left  with  a  problem  similar  to  the  multilevel  barium  cloud  striation 
problem  \  Lloyd  and  Haerendal,  1973;  Scannapieco  et  al. ,  1974;  1976;  Doles 
et  al.,  1976],  We  will  now  briefly  derive  the  multilevel  equations,  in 
general  form,  appropriate  to  our  ESF  problem. 

Consider  a  plasma  consisting  of  ions  and  electrons  imbedded  in  a  magnetic 
field  aligned  along  the  z  axis.  The  continuity  and  momentum  equations  de¬ 
scribing  the  system  are: 
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where  the  subscript  a  denotes  the  species  (i  for  ions,  e  for  electrons), 

n  is  the  species  number  density,  \/  is  velocity,  v  is  the  recombination  co- 

K 

efficient,  _E  is  the  electric  field,  £  is  the  gravitational  acceleration,  q 
is  the  species  charge,  v  is  the  species  collision  frequency  with  the  neutral 
atmosphere,  is  the  neutral  wind  velocity,  c  is  the  speed  of  light,  and 
m  is  the  species  mass.  Note  that  we  have  neglected  finite  temperature  effects 
(pressure),  and  the  effects  of  ion-ion  collisions  and  electron-ion  collisions 
(eventually,  we  will  even  neglect  electron-neutral  collisions).  We  further 
assume  that  we  are  interested  only  in  average  drift  velocities  over  time 
scales  long  compared  to  either  the  mean  time  between  collisions  or  the  gyro- 
period.  In  this  case  we  can  neglect  the  inertial  terms  (the  left  hand  side) 
of  (2),  and  invert  the  equation  to  obtain  an  algebraic  expression  for  v  : 
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The  vector  subscripts  x  and  jj refer  to  the  components  of  the  vector  which  are 

perpendicular  and  parallel  respectively  to  z.  We  take  ■  e  and  qg  -  -e. 

We  then  assume  that  v  /!J  *»  0  and  obtain 
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We  now  define  the  perpendicular  current 
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Substituting  (11)  through  (15),  (3)  and  (4)  into  (16), 

and  using  the  quasi- 

neutrality  approximation 
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For  our  problem 
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and  we  obtain 
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Since  0.01  iR^d.O  we  may  neglect  n^/R^  with  respect  to  m^. 

Defining  the  Pedersen  conductivity 
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and  noting  that  1-R^  1  =  -  v^/O^2  we  obtain 
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Equation  (23)  is~to~ oe  applied  to  each  of  our  layers  of  plasma.  Referring  to 
Figures  1  and  2  we  see  that  for  layers  1  and  3^=0  and  further  that 
_g  ■  -  Dgy  where  g  *  980  cm/sec  and  0<D<1  to  account  for  the  fact  that  £ 
is  not  perpendicular  to  15  for  plasma  away  from  the  equatorial  plane  (D  is 
actually  cos  6^  where  is  the  dip  angle).  Under  the  assumptions  we  shall 
make  later  it  will  be  seen  that  the  value  of  D  is  irrelevant  and  can  be  taken 


to  be  zero.  In  layer  2,  the  equatorial  plane,  we  have  £x  “  -  g  y,  and  we 

make  the  assumption  that  the  neutral  wind  is  directed  along  the  x  axis 

(U  ■  U  x,  where  x  points  west).  Furthermore,  since  layer  2  is  taken  to  be 
—  n 

at  F  region  altitudes  where  v^n/fl^<<l  (R^l ) ,  we  can  neglect  in  that  layer 
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the  second,  third,  and  fourth  terms  of  equation  (23)  with  respect  to  the  fifth, 
sixth,  and  first,  respectively.  So  we  have  for  the  three  layers: 


where  the  numerical  subscripts  refer  to  the  layer  numbers  depicted  in  Figure  2. 
Quasi-neutrality  demands  that 
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Integrating  (27)  along  field  lines  and  assuming  that  vanishes  at  z  =  +  “ 
we  obtain 
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If  we  model  our  plasma  as  an  array  of  discrete  layers  of  planes  of  plasma  per¬ 
pendicular  to  the  magnetic  field  as  in  Figure  2  we  may  replace  the  integral  i 

by  a  sum:  | 
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where  Az^  is  tfie  thickness  of  layer  k  measured  along  the  magnetic  field  line. 
By  our  assumption  of  equipotential  magnetic  field  lines  and  electrostatic 
electric  fields 


Eai  (x,y)  -  Ea2  (x,y)  -  E^  (x,y)  =  E  (x,y)  -  -  VJ  (x,y)  (31) 

where  we  have  neglected  the  slight  convergence  of  the  magnetic  field.  Then 
(30)  becomes 


where  the  subscript  b  denotes  the  sum  of  levels  1  and  3  and 
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(33) 


Also  we  have  defined 
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(34) 
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Note  Chat  implicit  in  the  above  manipulation  is  the  assumption  that  and 
are  constant  along  a  magnetic  field  line  within  a  given  layer,  as 
we  had  assumed  earlier.  Equations  (1)  and  (32)  constitute  the  system  of 
equations  we  must  solve.  In  general  it  will  be  necessary  to  resort  to 
numerical  means  for  this  task,  but  for  the  case  of  an  unperturbed  laminar 
ionosphere  it  is  both  possible  and  useful  to  find  a  simple  analytic  solution 
to  the  plasma  flow  field,  which  is  an  illuminating  example. 

Suppose  I  ,  E  ,  E  ,  are  functions  only  of  y  (altitude  in  the 
Pi  P2  P3 

equatorial  plane).  Then  for  any  set  of  boundary  conditions  on  <p  which  does 
not  itself  impose  an  x-dependence  on  <j>,  we  find  that  $  =  <Ky).  Then  (32) 
becomes 


i-  /< 
ay 


(I  +  E  +  E  )  (E  n.  —  u  ) 

pi  P2  P3  ay/  ay  P2  1  e  n' 


in. 


(35) 


the  general  solution  of  which  is 


in. 


(E  +E  +E  )  -  -  E_  -4  U_  +  J 


-  '  li  '  L. I  /  M 

pi  P2  P3  •  ay 


P2  i  e  n 


oy 


(36) 


where  is  a  constant,  and  we  have  dropped  the  subscript  2  on  m^, 


and 


3  $ 


Un.  Assuming  that  +  0  as  y  H  «  and  demanding  that  (or  equivalently 


P2 
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the  total  current)  vanish  at  y  «  +  “  we  get  -  0..  Recalling  that 

-  E  we  obtain 
9y  y 


y  E  +E  +E 

Pi  P2  P3 


a.  ^  "n 

i  e 


(37) 


The  ExB  plasma  motion  produced  by  this  electric  field  is  given  by 


cE 


v  = 

X 


E  +Z  +E  B 
Pi  P2  P3 


m. 

^  Q.  —  U 


(38) 


_E2- 


E  +E  +E 
Pi  P2  P3 


=  f  U 


where 


4 

c 


f  =  E  /  (E  +E  +E  ) 
P2  Pi  P2  P3 


(39) 


Note  that  the  plasma  drifts  at  a  fraction  f  of  the  neutral  wind  velocity,  and 
that  that  fraction  is  simply  the  ratio  of  the  "local"  (i.e.,  equatorial  plane) 
Pedersen  conductivity  to  the  total  field  line  conductivity  on  a  given  field 
line  (Note:  what  we  have  in  mind  here  and  in  the  numerical  simulations  is 
that  our  magnetic  field  line  integration  for  the  equatorial  F  region  is  over 
a  few  degrees  in  latitude,  and  that  regions  1  and  3  constitute  the  rest  of 
the  field  line  connected  ionosphere  as  a  load  on  the  circuit).  This  simple 
equation  has  some  remarkable  consequences  in  terms  of  the  motion  of  structures 
(spread  F  plumes,  for  example)  imbedded  in  the  equatorial  ionosphere.  Sup- 

IQdX 

pose  that  E  is  a  function  of  altitude  with  a  peak  E  at  altitude  h 

p2  P2  max 

Suppose  further  that  E  and  E  are  constants  such  that  E  +  E  *0.1  E1118*, 

Pi  P3  Pi  P3  P2 

and  that  we  impose  a  uniform  eastward  neutral  wind  of  100  m/sec  on  level  2 
(the  equatorial  plane).  We  now  create  a  model  ionosphere  (see  Table  1)  and 
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tabulate  the  eastward  plasma  velocity  as  a  function  of  altitude: 


TABLE  1 


Altitude  (km)  Z  /Z103* 

_  -22 _ E2_ 


Eastward  Plasma 
Velocity  (m/sec) 


600 

0.1 

50 

500 

0.5 

83 

400  (h  ) 
max 

1.0 

91 

300 

0.1 

50 

200 

0.01 

9 

Note  that  even  though  there  is  no  vertical  shear  in  the  neutral  wind  velocity, 

the  plasma  flow  field  contains  a  large  shear  with  opposing  signs  on  either 

side  of  h  .  The  effect  of  this  shear  is  to  bend  any  passive  vertical 
max 

structure  imbedded  in  this  flow  field  into  a  "C"  shape  as  depicted  in 

Figure  3  (Also  note  that  for  Z  *  Z  =0,  i.e.,  no  E  region,  from  (39) 

PI  P3 

f  *  1  and  the  plasma  moves  at  the  wind  speed  (Rishbeth,  1971)). 

The  above  result  is  quite  satisfying  in  that  it  offers  a  qualitative 
explanation  of  the  "C's",  "fishtails",  and  other  tilted  structures  seen  by 
Woodman  and  La  Hoa  [1976]  and  Tsunoda  [l98l]  in  their  observations  of  co¬ 
herent  radar  backscatter  from  the  meter-scale  irregularities  associated  with 
ESF  plumes.  However,  the  above  analysis  is  valid  only  for  passive  structures 
imbedded  in  a  laminar  unperturbed  ionosphere,  conditions  which  are  simply 
not  met  in  the  ESF  environment.  Numerical  simulations  are  necessary  to 
prove  the  case  unequivocally. 
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Fig.  3  —  Schematic  diagram  depicting  the  bending  of  a  passive  vertically-aligned 
structure  caught  up  in  a  velocity  shear  pattern  of  the  type  we  believe  exists  in  the 
nighttime  equatorial  ionosphere.  The  neutral  wind  is  eastward  and  uniform  in 
altitude,  and  the  response  of  the  plasma  (depicted  by  arrows)  is  to  move  at  some 
fraction  of  the  neutral  wind  velocity,  that  fraction  being  largest  at  the  altitude 
hmax  of  maximum  equatorial  plane  Pedersen  conductivity.  The  eastward  plasma 
velocity  falls  off  both  above  and  below  h_„,  as  shown. 


3.  Numerical  Simulations 

We  mentioned  before  that  equations  (1)  and  (32)  constitute  the  system 
we  wish  to  solve  numerically.  Let  us  now  be  more  specific.  Equation  (1) 
is  actually  six  equations  (one  electron  and  one  ion  equation  for  each  of  our 
three  layers).  By  quasi-neutrality  (17)  we  can  eliminate  three  of  these  and 
integrate  either  the  ion  or  electron  equation  at  each  level;  but  since  we 
have  made  the  assumption  that  =  0  (currents  along  field  lines  are  carried 
by  electrons)  we  can  more  conveniently  solve  the  two-dimensional  ion  continuity 
equation  at  each  layer: 

[if  +  7±  •  (n  vti)  -  -  vRn] 

k;  k  =  1,2,3  (40) 

In  the  simulations  we  present  here  we  have  set  vR  in  (40)  to  zero  for  simpli¬ 
city  and  because  at  the  F  region  altitudes  we  shall  be  dealing  with,  recombina¬ 
tion  effects  are  negligible. 

We  now  make  one  last  simplifying  assumption:  the  background  E  region 
plasma  (layers  1  and  3)  is  initially  uniform  in  density  and  Pedersen  conduc¬ 
tivity,  and  remains  so  during  the  course  of  our  simulation.  This  is  tantamount 
to  neglecting  compressibility  (Pedersen  mobility)  effects  in  the  E  region 
plasma.  Thus,  we  are  utilizing  layers  1  and  3  as  a  passive  load  in  an  iono¬ 
spheric  circuit,  in  order  to  allow  for  short  circuiting  effects.  A  true  multi¬ 
level  numerical  code  which  will  model  these  effects  self-consistantly  is  under 
development.  This  assumption  does  have  the  advantage,  though,  of  reducing  (40) 
to  a  single  equation  (since  3n/3t  -  0  for  levels  1  and  3)  and  of  eliminating 
H  and  the  last  two  terms  of  (32)  (since  all  the  terms  subscripted  by  b 
are  constant).  Our  final  pair  of  equations  to  be  solved  numerically  is  then 


(41) 


V  • 
x 


(n  ^1j? 


0 


(42) 


where  all  quantities  except  I  and  Z  refer  to  layer  2. 

Pi  P3 

Equation  (41)  is  solved  numerically  using  the  fully  multi-dimensional 
flux-corrected  transport  (FCT)  techniques  of  Zalesak  [1979].  Briefly, 

FCT  is  a  technique  originally  developed  by  Boris  and  Book  [l 97  3]  for  solving 
equations  of  the  form  (41)  where  steep  gradients  in  n  are  expected  to  form. 
The  fluxes  used  in  the  algorithm  are  nonlinear  weighted  averages  of  fluxes 
computed  by  high  and  low  order  finite  differences.  The  high  order  fluxes  are 
weighted  as  heavily  as  possible  subject  to  the  constraint  that  nonphysical 
oscillations  are  not  introduced.  Equation  (42)  is  solved  using  the  fully 
vectorized  incomplete  Cholesky  conjugate  gradient  (ICCG)  algorithm  of  Haln 
[l980],  which  is  an  extension  of  the  work  of  Kershaw  [l978].  This  algorithm 
is  extremely  fast  and  efficient  for  the  cases  described  below  for  which  the 
neutral  wind  was  set  to  zero;  however,  when  a  finite  eastward  neutral  wind 
was  used  the  ICCG  convergence  rate  became  painfully  slow  and  it  was  necessary 
to  resort  to  the  direct  elliptic  solver  of  Madala  [l 978]. 

The  numerical  calculations  to  be  presented  were  performed  on  a  two- 
dimensional  cartesian  mesh  using  40  points  in  the  x  (east-west)  direction  and 
140  points  in  the  y  (vertical  direction).  The  (uniform)  grid  spacing  was 
3  km  in  the  y  direction,  and  5  km  is  the  x  direction  for  all  calculations. 
Note  that  in  our  previous  work  we  used  2  km  spacing  in  the  y  direction. 

The  bottom  of  the  grid  corresponds  to  233  km  altitude  and  the  top  of  the 
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grid  to  676  km  altitude.  Periodic  boundary  conditions  were  imposed  on  both 
n  and  <f>  in  the  x  direction.  In  the  y  direction  transmissive  boundary  con¬ 
ditions  were  imposed  on  n  (3n/3x  =0),  and  the  normal  derivative  of  $  at  the 
top  boundary  was  chosen  such  that  the  normal  component  of  the  total  current 
(the  sum  over  all  three  layers)  was  zero  there  for  the  unperturbed  state. 

For  the  calculations  with  no  neutral  wind  this  implies  3<j>/3y  =0  at  the  upper 
boundary.  For  calculations  with  a  neutral  wind  this  implies 


f  +  [$«i  r  »„]  -  • 


(43) 


at  the  upper  boundary,  where  z°  is  the  Pedersen  conductivity  of  the  initial 
unperturbed  state.  At  the  lower  boundary  we  set  <p  =  0  for  all  cases. 

Three  kinds  of  plots  will  be  presented:  (1)  contours  of  constant 
n(x,y,t);  (2)  contours  of  constant  n(x,y,t)/nQ(y) ;  and  (3)  contours  of  con¬ 
stant  electrostatic  potential  <p(x,y,t).  Here  nQ(y)  is  the  initial  unperturbed 
electron  density  profile  in  layer  2.  Superimposed  on  each  contour  plot  is  a 
dashed  line  depicting  nQ(y)  for  reference  purposes.  Our  n^(y)  profile  is 
such  that  the  F£  peak  is  located  at  434  km  altitude,  and  the  minimum  electron 

density  scale  length  L  *  n  (3n  /3y)  1  is  10  km.  The  ion-neutral  collision 

o  o 

frequency  v^(y)  used  in  the  calculations  is  shown  in  Figure  4.  The  initial 
perturbation  used  to  start  each  calculation  was  a  mode  1  sine  wave  in  the 
x  direction: 


o^(y)  “  1  -  e  3  cos  (trx/100) 


(44) 


Three  calculations  were  performed  to  determine  the  effect  of  the 
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ALTITUDE  (km) 


0.01  0.1  I  10  100  200 


ION-NEUTRAL  COLLISION  FREQUENCY.  ujn(sec"‘) 

Fig.  4  —  The  ion-neutral  collision  frequency  (solid  line)  as  a  function 
of  altitude.  The  recombination  coefficient  was  set  to  zero  for  this 
study  (see  text). 
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background  E  region  plasma  and  of  the  eastward  neutral  wind  on  evolving 
spread  F  bubbles: 

1)  Calculation  2L,  in  which  £  =  £  =  U  =  0.  This  calculation  is 

Pi  P3  n 

identical  to  calculation  2L  of  Zalesak  and  Ossakow  [  198ol  except 
for  the  difference  in  vertical  grid  spacing  noted  previously. 

2)  Calculation  2LE,  identical  to  calculation  2L  above,  except  that  a 

constant  background  Pedersen  conductivity  has  been  included  such 
that  £  +  £  =  0.12  £°  ,  where  £°  is  the  maximum  Pedersen 

Pl  P3  2max  2 E03* 

conductivity  in  the  initial  unperturbed  equatorial  plane  (layer  2). 

We  believe  the  value  of  0.12  for  the  relative  background  Pedersen 
conductivity  level  to  be  a  conservative  figure.  Rishbeth  [l97l] 
used  a  value  of  0.2  as  being  representative  of  nighttime  conditions. 

3)  Calculation  2LEW,  identical  to  calculation  2LE  above,  except  that 
a  uniform  eastward  neutral  wind  of  150  m/sec  was  imposed  over  the 
entire  equatorial  plane  (U  =■  -  150  m/sec).  The  designations 

E  and  W  above  obviously  refer  to  the  presence  of  E  region  plasma  and 
neutral  winds,  respectively. 

Figure  5  shows  isodensity  contours  of  n(x,y,0)  for  our  initial  conditions 
(laminar  ionosphere  nQ(y)  plus  perturbation  (44)).  The  contours  are  labeled 
for  later  reference  purposes.  Note  that  in  this  and  all  subsequent  plots 
we  have  plotted  two  periods  (recall  that  we  have  periodic  boundary  conditions 
in  the  x  direction)  of  the  various  functions.  That  is,  the  40  by  140  mesh 
was  extended  to  80  by  140  for  plotting  purposes  only,  to  facilitate  comparison 
with  plots  of  calculations  run  with  a  neutral  wind,  in  which  structures  move 
across  the  grid. 

Figure  6  shows  isodensity  contours  of  n(x,y,t)  for  calculation  2L  at 
four  different  times  during  the  simulation.  Figure  7  shows  a  similar  sequence 
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ALTITUDE  (KM) 


EAST  (KM)  VEST 

Fig.  5  —  Iso-electron  density  contours  for  the  initial  perturbed  state  in  layer  2. 

This  represents  the  initial  conditions  for  our  numerical  simulation.  The  con¬ 
tours  are  labeled  in  units  of  electrons/cm3  in  E  format  notation  (1.0E1  =  1  X 
101 ,  etc.).  The  unperturbed  ionosphere  was  initially  laminar  (independent  of  x, 
the  east- west  direction)  and  is  exhibited  by  the  dashed  curve  showing  N0(y),  at 
any  point  in  the  east-west  (x)  direction.  This  curve  is  labeled  at  the  top  of  the 
figure.  The  perturbation  has  a  maximum  amplitude  of  e-3  in  relative  electron 
density,  is  a  pure  mode  1  sine  wave  in  x,  and  is  independent  of  altitude  y,  as 
described  in  the  text.  The  observer  is  looking  southward  so  that  B  is  out  of  the 
figure. 
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Fig.  7  —  Same  as  Figure  6,  but  for  calculation  2LE  (background  E  region, 
no  wind)  at  1500, 1900,  2215,  and  2560  seconds 


for  calculation  2LE.  (The  reader  may  note  that  the  2L  calculation  of  this 
paper  evolves  at  a  faster  rate  at  late  times  than  the  2L  calculation  of 
Zalesak  and  Os sakow  fl98ol.  This  is  primarily  due  to  two  improvements  in  the 
numerical  treatment  of  equations  (41)  and  (42),  implemented  since  our  previous 
studies:  1)  the  differencing  of  the  Hermitian  form  (42)  of  the  potential 
equation,  rather  than  the  non-Hermitian  expansion  we  were  constrained  to  use 
previously,  as  discussed*  in  the  appendix  of  Zalesak  and  Ossakow  [l98o3;  and 
2)  improved  treatment  of  the  continuity  equation  (41)  which  has  enabled  us  to 
further  reduce  the  numerical  diffusion  that  inevitably  occurs  across  electron 
density  gradients  as  steep  as  those  formed  at  the  edges  of  ESF  plumes  at  late 
times.  We  would  emphasize  that  the  conclusions  of  Zalesak  and  Ossakow  £l980} 
do  not  depend  on  late-time  rise  velocities  and  are,  therefore,  unaffected  by 
this  result).  In  comparing  Figure  6  with  Figure  7  we  are  looking  at  the  ef¬ 
fect  of  a  background  conductivity.  The  most  striking  difference  is  that  of 
the  time  scales,  whereby  2LE  takes  about  70Z  more  time  to  achieve  a  600  km 
altitude  plume  than  does  2L.  Qualitatively  this  can  be  understood  in  terms 
of  the  shorting  effect  of  the  background  E  region,  by  which  a  given  current 
can  be  driven  by  a  smaller  electric  field,  which  in  turn  means  smaller  plasma 
velocities.  Almost  as  striking  is  the  fact  that  the  2LE  plume  bifurcated 
while  the  2L  plume  did  not.  The  inevitability  of  the  bifurcation  in  calcula¬ 
tion  2LE  can  be  seen  even  in  the  very  early  time  plot  at  1500  sec,  where  the 
characteristic  flattening  of  a  significant  number  of  contours  in  the  upper 
portion  of  the  plume,  the  sure  signal  of  imminent  bifurcation  in  barium  cloud 
studies  T  Zabusky  et  al.*,  1973;  Scannapieco  et  al.,  1974;  Ossakow  et  al.,  1977; 
McDonald  et  al. ,  1980] ,  can  be  seen.  The  close  similarity  of  the  physics  of 
the  ESF  gravitational  instability  and  that  of  the  ExB  gradient  drift  in¬ 
stability  associated  with  the  bifurcation  and  striation  process  in  plasma 
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clouds,  has  been  noted  by  Scannapleco  and  Ossakow  [l976].  We  shall  draw 
heavily  on  our  knowledge  of  bifurcation  tendencies  in  plasma  clouds 
[*Ossakow  et  al . ,  1977;  McDonald  et  al.,  1980]  when  we  address  the  question  of 
why  the  plume  in  calculation  2LE  bifurcated  while  that  in  calculation  2L  did 
not,  later  in  this  paper. 

For  the  moment  we  note  that  there  are  two  primary  effects  of  the  pre¬ 
sence  of  a  background  conducting  region  (E  region):  1)  electric  fields  every¬ 
where  are  reduced  by  the  shorting  effects  of  the  background  conductivity,  re¬ 
sulting  in  an  overall  slower  evolution  for  the  instability;  and  2)  electric 
fields  are  reduced  the  most  in  the  regions  where  the  ratio  of  the  equatorial 
plane  conductivity  to  that  of  the  background  plane  is  smallest,  i.e.,  at  low 
altitudes,  rendering  the  2LE  configuration  incapable  of  drawing  plasma  from 
extremely  low  altitudes  to  produce  large  depletion  levels  inside  the  bubble. 
This  can  be  seen  easily  in  comparing  Figures  6  and  7  wherein  we  note  that 
the  isodensity  contours  at  low  altitudes  are  virtually  stationary  in  the  2LE 
case,  whereas  the  2L  configuration  results  in  significant  upward  movement  in 
even  the  lowest  density  plasma  near  the  bottom  of  the  plot.  The  more  effec¬ 
tive  shorting  of  the  electric  fields  by  the  background  layer  in  the  2LE 
case  is  seen  dramatically  in  Fig.  9a  and  9b,  where  we  plot  contours  of 
constant  electrostatic  potential  $  for  calculations  2L  and  2LE  respectively 
at  early  time.  (The  contour  level  increment  of  $  in  this  paper  is  chosen 
such  as  to  divide  the  maximum  excursion  of  <p  from  zero  into  7  equal  intervals 
The  contours  corresponding  to  positive  values  of  <|>  are  plotted  as  solid  lines 
while  those  corresponding  to  negative  values  are  plotted  as  dashed  lines. 

The  zero  contour  level  is  suppressed.  For  the  2L  and  2LE  cases  the 
symmetry  of  the  potential  would  cause  the  zero  contour  to  be  simply  two 
vertical  lines.  Since  the  electric  field,  and  hence  plasma  velocity,  is 
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inversely  proportional  to  the  contour  spacing,  this  normalization  allows  us 
to  easily  determine  by  eye  the  rate  at  which  the  upward  velocity  of  plasma 
is  decreasing  with  decreasing  altitude.  It  also  allows  us  to  visualize  the 
global  plasma  flow  field,  since  contours  of  $  are  essentially  streamlines  of 
the  plasma  flow).  Comparing  Figures  9a  and  9b,  we  note  a  much  more  rapid  de¬ 
crease  in  the  horizontal  component  of  the  electric  field  with  decreasing  al¬ 
titude  in  calculation  2L£  than  in  2L.  The  flow  field  in  2LE  is  mixing  plasma 
over  a  fairly  narrow  altitude  range,  while  that  in  2L  is  drawing  plasma  from 
deep  in  the  ionosphere,  where  the  plasma  densities  are  lowest.  Hence,  we 
should  expect  the  late  time  plume  in  calculation  2L  to  consist  of  plasma  of 
lower  density  (i.e.,  to  have  much  higher  depletion  levels)  than  that  in  cal¬ 
culation  2LE.  That  this  is  indeed  the  case  can  be  seen  in  Figures  10a  and 
10b,  where  we  compare  isodensity  contours  of  n(x,y)/nQ(y)  at  late  times  for 
the  two  calculations.  (Contours  of  n/nQ  in  this  paper  are  spaced  logarithmi¬ 
cally,  with  solid  lines  representing  depletions  (n/nQ<l)  and  dashed  lines  re¬ 
presenting  enhancements  (n/nQ>l).  The  kth  depletion  contour  represents  an 
n/nQ  value  of  0.5  ,  while  the  kth  enhancement  contour  represents  an  n/nQ 
value  of  2.0  ).  Although  the  bifurcation  of  the  2LE  plume  makes  the  compari¬ 
son  less  clean  than  it  would  otherwise  be,  it  is  obvious  that  the  depletion 
levels  of  the  2L  plume  are  much  higher  than  those  of  the  2LE  plume.  Depletion 
levels  (l-n/nQ)  in  the  upper  portions  of  the  2L  plume  are  greater  than  99. 2%, 
while  those  in  the  2LE  plume  are  only  about  94 1.  We  conclude  that  the  pre¬ 
sence  of  a  background  conducting  region  results  in  plumes  which  are  both 
slower  to  evolve  and  less  depleted  than  their  no-background  conterparts. 

In  Figure  8  we  present  isodensity  contours  of  calculation  2LEW  at  times 
similar  to  those  presented  for  the  2LE  calculation.  For  the  vin(y)  and  nQ(y) 
profiles  chosen  we  find  a  peak  in  E°  at  394  km,  40  km  below  the  F2  peak,  with 
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Fig.  9  —  Early  time  contours  of  constant  electrostatic  potential  0(x,y)  for  (a)  calculation 
2L  at  875  sec,  (b)  calculation  2  LE  at  1500  sec  and  (c)  calculation  2  LEW  at  1500  sec. 
The  potential  0o(y)  associated  with  the  unperturbed  initial  conditions  has  been  removed 
in  (c).  The  contour  level  increment  is  chosen  such  as  to  divide  the  maximum  excursion 
of  0  from  zero  into  7  equal  intervals.  The  contours  corresponding  to  positive  values  of  0 
are  plotted  as  solid  lines,  while  those  corresponding  to  negative  values  are  plotted  as  dashed 
lines.  The  zero  contour  level  is  suppressed. 


E°  falling  off  by  a  factor  of  ten  42  km  below  and  132  km  above  this  altitude 
P  2 

(at  352  km  and  526  km  altitude  respectively).  Using  Equation  (38)  to  approxi¬ 
mate  our  initial  shear  field  and  using  E  +  E  =  .12  I™3*  we  find  eastward 

Pi  P3  2 

plasma  drifts  of  134  m/sec  at  394  km  altitude  and  68  m/sec  at  both  352  and 
526  km  altitude.  The  plasma  shear  is  weaker,  but  over  a  larger  altitude  range, 
above  the  peak  in  than  below  it.  If  vertical  plasma  plumes  behave  as 
passive  structures,  we  would  expect  a  bending  of  the  structure  around  an 
altitude  of  394  km,  with  a  larger  slope  below  this  altitude  than  above  it. 
Looking  at  Figure  8  we  see  that  this  behavior  is  qualitatively  reproduced,  in 
spite  of  the  fact  that  the  self-consisteat  polarization  fields  produced  by 
the  plumes  represent  very  large  perturbations  on  the  equilibrium  fields  pro¬ 
ducing  the  plasma  shear.  In  Figure  8  we  have  placed  ourselves  in  a  frame 
moving  at  68  m/sec  eastward  to  minimize  both  computational  errors  and  computer 
time.  In  Figure  10c  we  show  late  time  isodensity  contours  of  n(x,y)/nQ(y) 
for  calculation  2LEW,  for  coraparision  to  Figures  IQa  and  10b.  The  bending 
of  the  plume  into  a  "C"  shape  about  an  altitude  of  360  km  is  quite  pro¬ 
nounced.  The  fact  that  this  "bending  point"  Is  more  than  30  km  below  the 
initial  maximum  in  equatorial  plane  Pedersen  conductivity  is  an  indication  of 
a  nonlinear  interaction  between  plume  rise  and  ambient  plasma  shear.  In  fact, 
this  shift  downward  can  be  understood  qualitatively  as  follows.  The  movement 
of  low  density  plasma  upward  inside  the  plume  is  accompanied  by  the  movement 
of  high  density  plasma  downward  in  the  regions  between  the  plumes  (see  Figures 
9  and  10) .  Since  the  scale  length  over  which  v  is  decreasing  with  altitude 
is  long  (^60  km)  compared  to  the  scale  length  over  which  the  electron  density 
is  increasing  with  altitude  (^10  km)  below  the  F2  peak,  the  effect  of  this 
downward  movement  of  high  density  plasma  is  to  move  the  point  of  maximum 
Pedersen  conductivity  in  the  equatorial  plane  and  hence  the  bending 
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Fig.  10  —  Late  time  contours  of  constant  n(x,y)/nG(y)  for  (a)  calculation  2L  at  1276  sec, 
(b)  calculation  2LE  at  2560  sec,  and  (c)  calculation  2LE  at  2331  sec.  The  contours  are 
spaced  logarithmically,  with  solid  lines  representing  depletions  (n/n„<  1)  and  dashed  lines 
representing  enhancements  (n/nc>  1).  The  kth  depletion  contour  represents  an  n/nQ  value 
of  0.5k,  while  the  kth  enhancement  contour  represents  an  n/n0  value  of  2.0k. 


point,  ‘downward. 

In  Figure  9  c,  we  show  contours  of  constant  <t> (x,y)-4>Q(y)  at  1500  sec 

for  calculation  2LEW,  for  comparison  to  Figures  9a  and  9b.  Here  4>0(y)  is 

the  initial  equilibrium  electrostatic  potential  of  the  unperturbed  initial 

conditions  (for  the  2L  and  2LE  cases  <f>o  =  0).  Subtracting  this  quantity 

from  <p  before  plotting  enables  us  to  examine  the  "underlying"  plume  motion 

on  which  the  shear  associated  with  the  initial  conditions  is  superposed. 

Remarkably,  the  motion  of  plasma  in  the  plume  is  not  purely  upward,  but 

rather  upward  and  westward,  despite  the  fact  that  we  have  removed  the 

asymmetry-inducing  profile  of  the  equilibrium  initial  wind  field.  Although 

this  simple  analysis  is  crude  (in  that  the  dependence  of  4>  on  the  plasma 

structure  is  not  linear,  i.e.,  in  Eq.  (42)  if  <f>A(x,y)  and  cf>B(x,y)  are 

solutions  for  EA  (x,y)  and  EB  (x,y)  respectively,  <(>A(x,y)  +  <j>B(x,y)  is  not 
P2  P2 

A  b 

a  solution  for  E  (x,y)  +  E  (x,y)),  it  would  seem  to  lend  support  to 
P2  P2 

the  ideas  advanced  by  Woodman  and  La  Hoz  [l976],  Ossakow  and  Chaturvedi  [1978], 
and  Ott  [l978]  who  proposed  that  a  neutral  wind  whose  eastward  velocity  ex¬ 
ceeded  that  of  the  plasma  would  combine  with  gravity  to  form  an  effective 
gravity  which  pointed  downward  and  eastward,  causing  bubbles  or  plumes  to 
drift  upward  and  westward  relative  to  the  surrounding  plasma.  Thus,  the 
westward  tilt  of  plumes  at  high  altitudes  would  appear  to  be  due  to  both  this 
effect  (since  we  have  shown  that  the  plasma  velocity  does  lag  the  neutral 
wind  velocity)  and  that  of  the  plasma  shear  which  we  have  addressed  earlier. 

We  would  point  out,  however,  that  the  mechanism  of  Woodman  and  La  Hoz  [l976], 
Ossakow  and  Chaturvedi  [l978],  and  Ott  [l978]  cannot  explain  the  eastward 
tilt  of  plumes  with  altitude  at  low  altitudes. 

In  comparing  calculations  2LE  and  2LEW  (Figures  7  and  8),  we  are 
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evaluating  the  effect  of  the  neutral  wind  itself,  since  the  same  background 
Pedersen  conductivity  was  used  in  both  cases.  We  wish  to  note  the  following 
points:  1)  the  two  calculations  evolve  at  approximately  the  same  rate  in 
time;  2)  the  primary  effect  of  the  wind  is  to  bend  the  plume  in  calculation 
2LEW  into  a  "C"  shape,  with  the  upper  part  of  the  "C"  being  much  larger  in 
altitude  extent  and  tilted  markedly  westward;  3)  the  plume  depletion  levels 
are  approximately  the  same  in  both  calculations;  and  4)  the  2LE  plume  bifur¬ 
cated  while  the  2 LEW  plume  did  not.  We  shall  address  this  last  question, 
along  with  the  question  of  why  the  plume  in  calculation  2L  did  not  bifurcate 
in  the  next  section. 
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4.  Discussion  and  Conclusions 

Before  proceeding  with  further  discussion  of  our  numerical  results,  let 
us  try  to  lend  support  to  the  idea  that  these  tilted  and  C-shaped  plumes 
are,  in  fact,  seen  in  equatorial  spread  F.  We  present  experimental  radar 
backscatter  maps  of  meter  scale  plasma  irregularities  taken  during  equatorial 
spread  F.  Figure  11  shows  a  map  of  3  meter  irregularities  (provided  by  J.P. 
McClure)  using  the  Jicamarca  radar.  Similar  plots  can  be  found  in  Woodman 
and  La  Hoz  f 1976]  .  Figure  12  shows  a  map  of  1  meter  irregularities  taken 
by  Tsunoda  [l98l]  using  the  ALTAIR  radar.  We  refer  the  reader  to  the  re¬ 
spective  papers  for  a  detailed  explanation  of  these  plots,  but  we  point  out 
that  the  Jicamarca  radar  scans  a  fixed  line  in  space,  and  plots  backscatter 
strength  as  a  function  of  time.  Therefore,  structures  caught  up  in  our 
postulated  plasma  shear  would  have  their  C-shaped  appearance  exaggerated  in 
the  Jicamarca  plots.  The  ALTAIR  radar,  however,  is  steerable;  and  it’s 
backscatter  plots  are  a  good  approximation  to  a  "snapshot"  of  the  backscatter 
strength  at  a  single  time.  In  both  plots  the  evidence  of  oppositely  tilting 
structures  at  high  and  low  altitudes  is  apparent.  In  making  comparisons  with 
these  small  scale  (^3m)  irregularity  radar  backscatter  maps  we  are  assuming 
that  these  maps  are  signatures  (Tsunoda,  1980)  of  the  large  scale  size  bubbles 
depicted  in  Figure  8.  That  is  the  steep  plasma  density  gradients  associated 
with  the  bubbles  in  Figure  8  drive  the  radar  backscatter  observed  irregulari¬ 
ties.  The  westward  and  upward  motion  of  the  bubbles  depicted  in  Figure  8 
is  in  agreement  with  the  satellite  in  situ  measurements  of  McClure  et  al., 
[1977]. 

It  is  our  belief  that  the  arguments  advanced  in  this  paper  offer  the 
most  plausible  explanation  yet  of  the  qualitative  behavior  of  equatorial 
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Fig.  12  ~  Maps  of  1  meter  irregularities  taken  by  the  steerable  ALT  AIR  radar  during 
equatorial  spread  F,  from  Tsunoda  [1981] .  These  maps  are  very  close  to  being  snap¬ 
shots  of  the  locus  of  irregularities  at  a  given  time.  Note  the  oppositely  tilting  struc¬ 
tures  at  different  altitudes  in  the  left-hand  map,  and  the  clear  bending  of  a  plume  into 
a  “C”  shape  about  an  altitude  of  about  400  km  in  the  right-hand  plot. 


spread  F  plumes:  "C"-shaped  structures  and  westward  plume  tilts  are  simply 
the  result  of  vertically  rising  spread  F  plumes  being  caught  up  in  the  ambient 
plasma  shear  both  as  they  rise  and  subsequently,  this  shear  being  the  natural 
consequence  of  a  neutral  wind  at  the  equator  and  E  regions  of  finite  conduc¬ 
tivity  connected  to  the  equatorial  F  region  along  magnetic  field  lines.  If 
there  is  a  neutral  wind,  but  no  E  region,  the  plasma  will  move  at  the  wind  speed, 
ESF  bubbles  will  rise  vertically,  and  the  attendant  radar  hackscatter  maps 
will  show  non-tilted  (i.e. ,  vertical)  plumes,  as  exhibited  in  the  measurements 
of  Kelley  et  al.,  (1981).  Furthermore,  even  without  the  equatorial  F  region 
neutral  wind  the  numerical  simulations  show  that  E  region  Pedersen  conductiv¬ 
ity  can  have  a  dramatic  effect  on  ESF.  For  example,  the  results  of  section 
3  show  that  ESF  has  been  slowed  down  as  has  the  attendant  bubble  evolution. 

In  addition,  the  ESF  bubble  in  the  presence  of  an  E  region  is  less  depleted 
than  without  the  E  region.  This  is  due  to  the  fact  that  the  E  region  has  a 
dramatic  effect  on  the  induced  polarization  electric  field  (see  Fig.  9)  which 
causes  the  rise  of  the  bubble,  the  fringe  field  component  of  which 
(Zalesak  and  Ossakow,  1980)  determines  the  region  below  the  F  peak  from 
which  plasma  is  drawn  (i.e.,  which  plasma  makes  up  the  bubble).  The  pre¬ 
sence  or  absence  of  an  E  region  could  also  explain  why  bubbles  (with  large 
depletions)  stop  rising  at  altitudes  of  400-500  km  (see  McClure  et  al.,  1977) 
on  some  occasions,  but  not  on  others.  Indeed,  in  addition  to  the  height  of 
the  F  peak  and  bottomside  background  electron  density  gradient  scale  lengths 
(Ossakow  et  al.,  1979)  influencing  ESF  evolution,  the  E  region  conductivity 
could  also  determine  why  one  does  or  does  not  observe  ESF  even  when  the 
previously  mentioned  conditions  are  satisfied.  Thus,  the  E  region  (even  at 
night)  could  be  a  controlling  factor  in  the  formation  of  ESF  irregularities 


and  could  account  for  such  things  as  the  longitudinal  influence  (JBasu  and 
Kelley,  1977;  Livingston,  1980)  on  ESF  formation  and  phenomena. 

There  are  two  matters,  however,  which  bear  further  discussion:  1)  the 
question  of  why  the  plume  in  calculation  2LE  bifurcated  while  those  in  cal¬ 
culations  2L  and  2LEW  did  not;  and  (2)  the  question  of  where  along  the  edge  of 
the  primary  plumes  in  calculation  2LEW  one  should  expect  to  see  secondary 
instabilities.  If  we  note  that  an  equatorial  spread  F  plume  (bubble)  is 
nothing  more  than  an  inverse  plasma  cloud  "finger"  (i.e.,  an  elongated  region 
of  low  density  plasma  (an  ESF  "plume")  penetrating  a  region  of  high  density 
plasma  is  the  inverse  of  an  elongated  region  of  high  density  plasma  (a  plasma 
cloud  "finger")  penetrating  a  region  of  low  density  plasma),  we  find  that 
the  question  of  why  the  2LE  plume  bifurcated  and  the  2L  plume  did  not  has 
already  been  answered  for  us.  McDonald  et  al.,  [l 98 l]  have  shown  in  their 
study  of  bifurcation  tendencies  of  plasma  cloud  fingers  that  the  critical 
quantity  determining  the  speed  with  which  a  plasma  finger  will  bifurcate  is 
M,  the  ratio  of  the  Pedersen  conductivities  inside  and  outside  the  finger. 
When  M  is  moderate,  in  the  range  2  to  10,  the  bifurcation  tendency  is  high, 
while  for  M  near  1  or  M  greater  than  100,  the  bifurcation  tendency  is  ex¬ 
tremely  small.  Looking  at  the  2LE  plume  (Fig.  7),  and  recalling  that  we 
have  a  background  Pedersen  conductivity  of  0.12  times  the  maximum  equatorial 
plane  Pedersen  conductivity  we  find  that  M  *  (the  relevant  quantity  since  we 
are  dealing  with  inverse  plasma  clouds  here)  in  the  2LE  plume  is  about  9, 
making  it  a  prime  candidate  for  bifurcation;  while  M-1  for  the  2L  plume 
(Fig.  6)  is  about  104,  indicating  a  bifurcation  tendency  near  zero.  The 
question  of  why  the  2LEW  plume  did  not  bifurcate  is  a  little  harder  to 
answer.  Based  on  the  arguments  advanced  above,  the  2LEW  plume  would  have 
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been  just  as  likely  a  candidate  for  bifurcation  as  the  2LE  plume.  However, 
we  note  that  the  2LEW  plume  is  rising  into  a  region  of  very  strong  plasma 
shear.  Perkins  and  Doles  Ql 9 7 5^  have  shown  that  such  a  shear  would  provide 
a  stabilizing  mechanism  for  any  secondary  instability  (i.e.,  bifurcation) 
which  attempted  to  grow  on  the  topside  of  the  2 LEW  plume,  although  the  geome¬ 
try  used  in  their  study  was  considerably  simpler  than  that  associated  with 
a  rising  ESF  plume.  We  advance  this  mechanism  as  a  plausible,  but  less  than 
totally  convincing, explanation  of  why  the  2LEW  plume  did  not  bifurcate,  only 
because  we  can  provide  no  other  at  this  time. 

The  question  of  secondary  instabilities  on  the  perimeters  or  in  the 
interiors  of  the  plumes  is  in  many  ways  the  most  interesting  aspect  of  this 
study.  We  will  confine  ourselves  to  the  2LEW  plume  (Fig.  8)  at  late  time 
(2331  sec)  since  it  is  both  the  most  interesting  and  the  most  realistic. 

Given  the  limited  spatial  resolution  of  these  numerical  studies,  it  must  be 
realized  that  the  actual  numerical  simulation  of  the  evolution  of  small-scale 
secondary  instabilities,  within  the  context  of  the  present  simulations,  is 
an  impossibility.  However,  we  do  have  the  resolution  to  be  able  to  observe 
the  precursor  of  the  plasma  fluid  instabilities  that  we  believe  are  active: 
the  steepening  of  electron  density  gradients.  By  observing  the  location 
of  these  regions  of  steepening,  and  by  augmenting  this  procedure  with  a  cell 
by  cell  local  stability  analysis,  we  should  be  able  to  predict  both  the 
location  of  secondary  instabilities  and  the  mechanism  causing  them.  We  use 
the  term  "local  stability  analysis"  to  mean  a  local  evaluation,  numerical  in 
this  case,  of  the  generalized  gradient  drift  growth  rate  y  given  by 
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Note  that  E  in  Eq.  (45)  above  includes  the  self-consistant  polarization 
electric  field  given  by  the  solution  to  the  potential  equation  (42).  The 
influence  of  this  term  on  y  is  large,  and  any  stability  analysis  which  were 
to  ignore  these  polarization  fields  would  rest  on  shaky  ground.  Among 
the  assumptions  implicit  in  the  application  of  Eq.  (45)  are:  1)  the  growth 
rates  y  are  large  compared  to  the  speed  with  which  the  primary  Rayleigh- 
Taylor  mode  is  evolving;  and  2)  the  k-vector  associated  with  the  growing 
perturbation  is  perpendicular  to  (which  gives  maximum  growth). 

Looking  at  Figure  8  at  the  latest  time  (2331  sec)  we  see  that  the  pri¬ 
mary  regions  of  steepening  are  two:  the  west  wall  of  the  plume  at  low  al¬ 
titudes  (below  370  km  for  this  particular  model)  and  the  east  wall  of  the 
plume  at  higher  altitudes  (above  370  km).  Local  stability  analysis  verifies 
that  these  are  precisely  the  regions  of  largest  growth  rate  for  the  gradient- 
drift/gravitational  Rayleigh-Taylor  instability.  A  less  complete  analysis  of 
just  the  effect  of  an  eastward  neutral  wind  on  a  more  or  less  vertical  plume 
would  predict  that  only  the  west  wall  of  the  plume  would  be  unstable,  but 
this  analysis  neglects  the  effects  of  the  bending  of  the  plume  which  orients 
the  normally  stable  east  wall  of  the  plume  so  that  it  is  once  again  unstable 
to  the  gravitational  instability,  and  the  effects  of  the  polarization  electric 
field  produced  self-consistently  by  the  ionosphere-plume  system,  whose  effect 
is  to  mitigate  the  expected  wind-driven  instability  over  most  of  the  plume 
structure.  We  conclude,  then,  that  for  this  particular  plume,  we  would  ex¬ 
pect  secondary  instabilities  along  the  west  wall  at  low  altitudes  and  along 
the  east  wall  at  high  altitudes,  with  the  "switch"  taking  place  at  about 
370  km  altitude.  If  these  instabilities  eventually  cascade  down  to  smaller 
and  smaller  scale  sizes  (or  provide  the  steep  plasma  density  gradients 
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necessary  for  further  instability  mechanisms),  eventually  reaching  the  1-3 
meter  scale  sizes  seen  on  backscatter  radar  maps,  we  would  expect  the  radar 
maps  to  trace  out  the  locus  of  the  west  wall  of  the  plume  at  low  altitudes, 
and  the  east  wall  at  high  altitudes,  giving  rise  to  an  even  more  exaggerated 
"C"  shape  than  that  of  the  simulation  plume  (bubble)  itself  (see  the  ex¬ 
aggerated  "C"  traced  out  by  the  locus  of  steepened  gradients  in  Figure  8 
at  2331  sec). 

We  wish  to  close  this  section  by  briefly  reviewing  work  by  other  re¬ 
searchers  which  we  believe  has  relevance  to  the  results  presented  here. 

Two  recent  papers  have  shown  experimental  evidence  of  a  shear  in  east-west 
plasma  motion  in  the  equatorial  ionosphere:  Kudeki  et  al . ,  [  198l3  and 
Tsunoda  et  al. ,  [l98]]  .  Both  papers  show  evidence  of  an  increase  in  eastward 
plasma  velocity  with  altitude,  in  agreement  with  the  behavior  we  postulate 
here  for  altitudes  below  the  peak  in  equatorial  F  region  integrated  Pedersen 
conductivity.  It  is  our  belief  that  experimental  observations  at  even  higher 
altitudes  than  that  examined  in  the  above  papers  would  show  a  decrease  in 
eastward  plasma  velocity  with  altitude  at  these  higher  altitudes  (and 
possibly  even  westward  velocities),  In  a  manner  similar  to  that  shown  in 
Figure  3.  Both  Kudeki  et  al . ,  [l98l]  and  Tsunoda  et  al.,[l98l]  show 
evidence  of  a  plasma  velocity  reversal  point,  that  is,  an  altitude  below  the 
F2  peak  below  which  the  plasma  velocity  actually  becomes  westward  (as  the 
eastward  velocity  passes  through  zero).  The  simple  model  we  have  presented 
here  offers  no  explanation  for  this  phenomenon.  The  reason  is  that  we  have 
assumed  here  that  the  E  regions  connected  to  the  equatorial  F  region  along 
field  lines  are  passive  and  free  of  any  dynamics  of  their  own.  Actually, 
however,  these  E  regions  are  subject  to  strong  diurnal  tidal  neutral  winds, 
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which  are  westward  at  the  times  associated  with  spread  F.  At  equatorial  al¬ 
titudes  well  below  the  F2  peak,  the  field-line  integrated  Pedersen  conductiv¬ 
ity  is  dominated  by  that  of  the  E  regions,  and  hence  the  westward  neutral 
winds  in  the  E  regions  are  able  to  set  up  polarization  electric  fields 
which  impress  themselves  on  the  equatorial  region,  causing  a  corresponding 
westward  drift  of  plasma  in  the  equatorial  plane  at  low  altitudes. 

In  fact,  had  we  retained  the  neutral  wind  terms  in  layers  1  and  3,  Eq. 
(38)  would  have  become 
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where  U  ,  U  ,  and  U  are  the  east-west  neutral  winds  in  layers  1,  2, 
ni  n2  n3 

and  3  respectively.  If  U  is  eastward  and  both  U  and  U  are  westward, 

n2  ni  n3 

it  is  obvious  that  westward  plasma  velocities  will  exist  at  any  altitude 
for  which 

Z  U  4-  Z  U  , 

.  £3.  .S3 _ El_ni  (a 

U  k 

n2 

This  effect  and  the  consequent  plasma  velocity  reversal  point  were  first 
described  by  Hellis  et  al. , [1974],  whose  detailed  self-consistent  numerical 
model  of  the  E-and  F-region  neutral  gas  and  plasma  system  also  shows  both 
the  plasma  shear  and  an  altitude  at  which  the  eastward  plasma  velocity 
maximizes,  as  we  have  proposed  here.  In  fact,  at  extreme  equatorial  F 
region  altitudes,  plasma  well  away  from  the  equatorial  plane  (both  E  and  F 
region  plasma)  may  again  dominate  the  integrated  Pedersen  conductivity,  and 
if  the  corresponding  neutral  winds  are  westward,  we  should  expect  to  see 
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westward  plasma  drifts  at  these  high  altitudes,  as  mentioned  above. 

The  influence  of  contributions  to  field  line  integrated  conductivity 
from  plasma  away  from  the  equatorial  plane  on  the  rise  of  ESF  bubbles  has 
also  been  addressed  by  Anderson  and  Haerendel  [l979].  In  their  model  they 
incorporated  flux  tube  integrated  quantities  of  electron  content  and  Pedersen 
conductivity  and  utilized  a  one-dimensional  sheet  model  for  the  bubble. 

This  latter  assumption  resulted  in  a  simple  algebraic  expression  for  the 
induced  polarization  electric  field  inside  the  bubble  in  terms  of  the  flux 
tube  integrated  quantities.  ESF  bubble  rise  in  the  collision  dominated 
Rayleigh- Taylor  regime  with  and  without  an  ambient  eastward  electric  field, 

E^,  was  investigated.  These  authors  noted  that  flux  tube  integrated  quanti¬ 
ties  could  have  a  considerable  influence  on  the  outcome  of  ESF  bubble  rise,  an 
observation  consistent  with  our  comparison  of  the  2L  and  2LE  cases.  Burke 
[l979]  analyzed  the  effect  of  the  sunrise  "turning-on"  of  the  E  region  and  its 
subsequent  contribution  to  the  demise  of  ESF  bubbles  in  the  topside  F  region 
ionosphere  near  the  dawn  terminator.  A  simple  analytic  model  was  used  which 
showed  how  electric  fields  within  the  bubbles  could  be  discharged  through 
the  conducting  sunrise  E  region.  The  results  showed  that  the  conducting 
E  region  could  effectively  halt  the  upward  bubble  rise  velocity.  This  im¬ 
portance  of  an  E  region  is  consistent  with  our  findings  in  comparing  the  2L 


and  2LE  cases. 


** 


5.  Future  Work 


The  "three  layer  model"  used  in  this  study  is  simple  to  be  sure,  al¬ 
though  we  believe  it  adequately  describes  the  qualitative  behaviour  of  ESP 
plumes  in  terms  of  C-shaped  structures,  westward  tilts,  and  the  effects  of 
a  background  E  region .  Work  toward  improving  the  model  and  its  input  is 

ongoing  on  several  fronts.  First,  we  would  like  to  make  the  E  regions 
"active",  i.e.,  to  allow  external  forces,  such  as  neutral  winds  to  act  on 
the  E  region  plasma,  and  to  self-consistently  solve  the  continuity  equation 
there.  Second,  we  would  like  to  better  resolve  the  plasma  distribution  along 
magnetic  field  lines  by  adding  more  layers  to  the  model  to  represent  plasma 
between  the  equatorial  plane  and  the  E  regions.  The  total  number  of  layers 
might  be  seven  or  nine.  Third,  we  would  like  to  incorporate  more  realistic 
models  of  electron  and  neutral  density  distributions,  external  electric  fields, 
neutral  winds,  and  chemistry  into  the  models.  Sources  for  this  information 
might  be  empirical  data  or  models  of  the  type  developed  by  Anderson  [1973 ] 
and  by  Forbes  and  Garrett  [l978] .  Last,  but  certainly  not  least,  a  continuing 
effort  is  being  made  to  keep  the  numerical  techniques  used  in  the  code  as 
close  to  state-of-the-art  as  possible.  A  recent  advance  [zalesak.  198 1] 
should  significantly  improve  our  already  quite  good,  but  certainly  not 
perfect,  numerical  algorithms  for  solving  the  continuity  equation  in  the 


very  near  future. 
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Dr.  H.  B.  HcElroy 
Dr.  R.  Llndzen 

Pennsylvania  State  University 
University  Park,  Pennsylvania  16802 

Or.  J.  S.  Nlsbet 
Or.  P.  R.  Rohrbeugh 
Dr.  L.  A.  Carpenter 
Dr.  H.  Lee 
Or.  R.  Olvany 
Or.  P.  Bennett 
Or.  E.  Kievans 


University  of  California,  Berkeley 
Berkeley,  California  94720 

Or.  M.  Hudson 

Utah  State  University 
4th  N.  and  8th  Streets 
Logan,  Utah  84322 

Or.  P.  H.  Banks 
Or.  R.  Harris 
Dr.  K.  Baker 

Cornell  University 
Ithaca,  New  fork  14850 

Dr.  W.  E.  Swartz 
Or.  R.  Sudan 
Dr.  0.  Farley 
Dr.  N.  Kelley 


NASA 

Goddard  Space  Flight  Center 
Greenbelt,  Nary  land  20771 

Or.  S.  Chandra 
Or.  K.  Maede 
Dr.  R.F.  Benson 

Princeton  University 
Plasma  Physics  Laboratory 
Princeton,  New  Jersey  08540 

Dr.  F.  Perkins 
Or.  E.  Frlemin 

Institute  for  Defense  Analysis 
400  Army/Navy  Drive 
Arlington,  Virginia  22202 

Dr.  E.  Bauer 

University  of  Narytand 
College  Park,  HD  20742 

Dr.  K.  Papadopoulos 
Dr.  E.  Ott 

University  of  Pittsburgh 
Pittsburgh,  Pa.  l^>2ol 

Or.  N.  Zabusky 
Dr.  H.  Blondl 

Defense  Documentation  Center 
Cameron  Station 
Alexandria.  VA  22314 

(12  copies  If  open  publication 
otherwise  2  copies) 

12CY  Attn  TC 

University  of  California 

Los  Alamos  Scientific  Laboratory 

J-10,  NS- 664 

Los  Alamos,  New  Hexlco  87545 

N.  Ponqracz 
D.  Simons 
G.  Baresch 
L.  Duncan 

Massachusetts  Institute  of  Technology 
Plasma  Fusion  Center 
Library,  NW16-262 
Cambridge.  HA  02139 


University  of  California.  Los  Angeles 
105  Hlllgard  Avenue 
os  Angeles,  California  90024 

Or.  F.  V.  Coronltl 
Or.  C.  Kennel 


57 


